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Abstract 

Determinant computation is the core procedure in many important ge- 
ometric algorithms, such as convex hull computations and point locations. 
As the dimension of the computation space grows, a higher percentage of 
the computation time is consumed by these predicates. In this paper we 
study the sequences of determinants that appear in geometric algorithms. 
We use dynamic determinant algorithms to speed-up the computation of 
each predicate by using information from previously computed predicates. 

We propose two dynamic determinant algorithms with quadratic com- 
plexity when employed in convex hull computations, and with linear com- 
plexity when used in point location problems. Moreover, we implement 
them and perform an experimental analysis. Our implementations out- 
perform the state-of-the-art determinant and convex hull implementations 
in most of the tested scenarios, as well as giving a speed-up of 78 times 
in point location problems. 

Keywords: computational geometry, determinant algorithms, orienta- 
tion predicate, convex hull, point location, experimental analysis 

1 Introduction 

Detcrminantal predicates are in the core of many important geometric algo- 
rithms. Convex hull and regular triangulation algorithms use orientation pred- 
icates, the Delaunay triangulation algorithms also involve the in-sphere predi- 
cate. Moreover, algorithms for exact volume computation of a convex polytope 
rely on determinantal volume formulas. In general dimension d, the orientation 
predicate of d+ 1 points is the sign of the determinant of a matrix containing the 
homogeneous coordinates of the points as columns. In a similar way, the volume 
determinant formula and in-sphere predicate of d + 1 and d + 2 points respec- 
tively can be denned. In practice, as the dimension grows, a higher percentage 
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of the computation time is consumed by these core procedures. For this reason, 
we focus on algorithms and implementations for the exact computation of the 
determinant. We give particular emphasis to division-free algorithms. Avoiding 
divisions is crucial when working on a ring that is not a field, e.g., integers 
or polynomials. Determinants of matrices whose elements are in a ring arise 
in combinatorial problems [Kra05 , in algorithms for lattice polyhedra [BP99 
and secondary polytopes iRam02| or in computational algebraic geometry prob- 
lems |CLO05| . 

Our main observation is that, in a sequence of computations of determinants 
that appear in geometric algorithms, the computation of one predicate can be 
accelerated by using information from the computation of previously computed 
predicates. In this paper, we study orientation predicates that appear in convex 
hull computations. The convex hull problem is probably the most fundamental 
problem in discrete computational geometry. In fact, the problems of regular 
and Delaunay triangulations reduce to it. 

Our main contribution is twofold. First, we propose an algorithm with qua- 
dratic complexity for the determinants involved in a convex hull computation 
and linear complexity for those involved in point location problems. Moreover, 
we nominate a variant of this algorithm that can perform computations over the 
integers. Second, we implement our proposed algorithms along with division- 
free determinant algorithms from the literature. We perform an experimental 
analysis of the current state-of-the-art packages for exact determinant computa- 
tions along with our implementations. Without taking the dynamic algorithms 
into account, the experiments serve as a case study of the best implementation 
of determinant algorithms, which is of independent interest. However, dynamic 
algorithms outperform the other determinant implementations in almost all the 
cases. Moreover, we implement our method on top of the convex hull package 
triangulation BDH09J and experimentally show that it attains a speed-up up 
to 3.5 times, results in a faster than state-of-the-art convex hull package and a 
competitive implementation for exact volume computation, as well as giving a 
speed-up of 78 times in point location problems. 

Let us review previous work. There is a variety of algorithms and imple- 
mentations for computing the determinant of a d x d matrix. By denoting 
0(d u ') their complexity, the best current ui is 2.697263 |KV05j . However, good 
asymptotic complexity does not imply good behavior in practice for small and 
medium dimensions. For instance, LinBox [DGG + 02] which implements al- 
gorithms with state-of-the-art asymptotic complexity, introduces a significant 
overhead in medium dimensions, and seems most suitable in very high dimen- 
sions (typically > 100). Eigen [GJ + 10 and CGAL |CGA] implement decompo- 
sition methods of complexity 0(n 3 ) and seem to be suitable for low to medium 
dimensions. There exist algorithms that avoid divisions such as |Rot01j with 
complexity 0(n i ) and [Birll with complexity 0{nM{n)) where M(n) is the 
complexity of matrix multiplication. In addition, there exists a variety of algo- 
rithms for determinant sign computation BEPP99, ABM99] . The problem of 
computation of several determinants has also been studied. TOPCOM |Ram02] . 
the reference software for computing triangulations of a set of points, efficiently 
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precomputes all orientation determinants that will be needed in the computa- 
tion and stores their signs. In [EFKP12 , a similar problem is studied in the 
context of computational algebraic geometry. The computation of orientation 
predicates is accelerated by maintaining a hash table of computed minors of the 
determinants. These minors appear many times in the computation. Although, 
applying that method to the convex hull computation does not lead to a more 
efficient algorithm. 

Our main tools are the Sherman-Morrison formulas [SM50, Bar51j. They 
relate the inverse of a matrix after a small-rank perturbation to the inverse of 
the original matrix. In |San04j these formulas are used in a similar way to solve 
the dynamic transitive closure problem in graphs. 

The paper is organized as follows. Sect. [2] introduces the dynamic deter- 
minant algorithms and the following section presents their application to the 
convex hull problem. Sect. [4] discusses the implementation, experiments, and 
comparison with other software. We conclude with future work. 



2 Dynamic Determinant Computations 

In the dynamic determinant problem, a d x d matrix A is given. Allowing 
some preprocessing, we should be able to handle updates of elements of A and 
return the current value of the determinant. We consider here only non-singular 
updates, i.e., updates that do not make A singular. Let (A)i denote the i-th 
column of A, and e% the vector with 1 in its i-th place and everywhere else. 

Consider the matrix A' , resulting from replacing the i-th column of A by 
a vector u. The Sherman-Morrison formula SM50, B ar5T] states that (A + 

wv T )~ 1 = A~ x — ^jt^a=t^ ^ ■ An i-th column update of A is performed by 
substituting v — e% and w — u — (A)i in the above formula. Then, we can write 
A' -1 as follows. 



If A 1 is computed, we compute A 1 1 using Eq.[l]in 3d 2 + 2d + 0(l) arithmetic 
operations. Similarly, the matrix determinant lemma |Har97] gives Eq. [2] below 
to compute det(A') in 2d + 0(1) arithmetic operations, if det(A) is computed. 

det(A') = det{A + (u- {A)i)ef) = (1 + ef A-\u - (A),)det{A) (2) 

Eqs.[l]and[2]lead to the following result. 

Proposition 1 |SM50j The dyn amic determinant problem can be solved using 
0{d u) ) arithmetic operations for preprocessing and 0(d 2 ) for non-singular one 
column updates. 

Indeed, this computation can also be performed over a ring. To this end, 
we use the adjoint of A, denoted by A adj , rather than the inverse. It holds that 
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A ad i — det(A)A 1 , thus we obtain the following two equations. 



A' adj = fA A A adj det(A') - (A adj (u - (A),)) [ej A adj )) (3) 

det(A') = det(A) + ej A adj (u ~ (A) t ) (4) 

The only division, in Eq.[3J is known to be exact, i.e., its remainder is zero. The 
above computations can be performed in 5c? 2 + d + O(l) arithmetic operations 
for Eq. [3] and in 2c? + 0(1) for Eq. |4j In the sequel, we will call dynAnv the 
dynamic determinant algorithm which uses Eqs. [T] and [2] and dyn_adj the one 
which uses Eqs. [3] and |4j 



3 Geometric Algorithms 

We introduce in this section our methods for optimizing the computation of se- 
quences of determinants that appear in geometric algorithms. First, we use dy- 
namic determinant computations in incremental convex hull algorithms. Then, 
we show how this solution can be extended to point location in triangulations. 

Let us start with some basic definitions from discrete and computational 
geometry. Let A C K d be a pointset. We define the convex hull of a pointsct 
A, denoted by conv(„4), as the smallest convex set containing A. A hyper- 
plane supports conv(A) if conv(.A) is entirely contained in one of the two closed 
half-spaces determined by the hyperplane and has at least one point on the 
hyperplane. A face of conv(A) is the intersection of conv(.4) with a support- 
ing hyperplane which does not contain conv(„4). Faces of dimension 0, l,d — 1 
are called vertices, edges and facets respectively. We call a face / of conv(^4) 
visible from a S R d if there is a supporting hyperplane of / such that conv(^4) 
is contained in one of the two closed half-spaces determined by the hyperplane 
and a in the other. A fc-simplex of A is an affinely independent subset S of 
A, where dim(conv(S)) = k. A triangulation of A is a collection of subsets of 
A, the cells of the triangulation, such that the union of the cells' convex hulls 
equals conv(„4), every pair of convex hulls of cells intersect at a common face 
and every cell is a simplex. 

Denote a the vector (a, 1) for a G R d . For any sequence C of points € 
A, i = 1 . . . d + 1, we denote Ac its orientation (d + 1) x (d + 1) matrix. For 
every a^, the column i of Ac contains a^'s coordinates as entries. For simplicity, 
we assume general position of A and focus on the Beneath-and-Beyond (BB) 
algorithm Sci81 . However, our method can be extended to handle degenerate 
inputs as in [Ede871 Sect. 8.4], as well as to be applied to any incremental 
convex hull algorithm by utilizing the dynamic determinant computations to 
answer the predicates appearing in point location (see Cor.[2|. In what follows, 
we use the dynamic determinant algorithm dyn^adj, which can be replaced by 
dyn_inv yielding a variant of the presented convex hull algorithm. 

The BB algorithm is initialized by computing a c?-simplex of A. At every 
subsequent step, a new point from A is inserted, while keeping a triangulated 
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Algorithm 1: Incremental Convex Hull (A) 

Input : pointset A C R d 
Output: convex hull of A 

sort A by increasing lexicographic order of coordinates, i.e., 
A = {ax, . . . ,a n }; 

T { d-face of conv(a\, a^+i)}; // conv{a\, . . . , a^+i) is a 

rf-simplex 
Q <— { facets of conv(a\, . . . , cid+i)}; 
compute Af* u ... >ad+l} ,det(A {au ^ ad+l} ); 
foreach a e {a d+2 , ... , a n } do 

Q'^Q; 

foreach F £ Q do 

C <— the unique d-face s.t. CeT and F e C; 

u <— the unique vertex s.t. wfC and u ^ F; 

C'^F\j{a}; 

i <— the index of u in Ac; 

// both det(Ac) and AJv 7 were computed in a previous 
step 

det(Ac) det(Ac) + {A a * J ) l {u - a); 

if det(A C /)det(A c ) < and &et(A c >) ^ then 

T<-TU {d-face of conu(C')}; 

Q' ■<— Q' G {(d — l)-faces of C"}; // symmetric difference 

_Q^Q'; 
return Q; 



convex hull of the inserted points. Let t be the number of cells of this trian- 
gulation. Assume that, at some step, a new point a £ A is inserted and T is 
the triangulation of the convex hull of the points of A inserted up to now. To 
determine if a facet F is visible from a, an orientation predicate involving a and 
the points of F has to be computed. This can be done by using Eq.|4]if we know 
the adjoint matrix of points of the cell that contains F. But, if F is visible, 
this cell is unique and we can map it to the adjoint matrix corresponding to its 
points. 

Our method (Alg. [T]), as initialization, computes from scratch the adjoint 
matrix that corresponds to the initial d-simplex. At every incremental step, 
it computes the orientation predicates using the adjoint matrices computed in 
previous steps and Eq. [4] It also computes the adjoint matrices corresponding 
to the new cells using Eq. [3j By Prop. [T] this method leads to the following 
result. 

Theorem 1 Given a d- dimensional pointset all, except the first, orientation 
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predicates of incremental convex hull algorithms can be computed in 0(d 2 ) time 
andO(d 2 t) space, where t is the number of cells of the constructed triangulation. 

Essentially, this result improves the computational complexity of the deter- 
minants involved in incremental convex hull algorithms from 0(d UJ ) to 0(d 2 ). 
To analyze the complexity of Alg. [Tl we bound the number of facets of Q in 
every step of the outer loop of AlgTTT] with the number of (d — l)-faces of the 
constructed triangulation of conv(_4), which is bounded by (d + l)t. Thus, using 
Thm. [T] we have the following complexity bound for Alg. [T] 

Corollary 1 Given n d-dimensional points, the complexity of BB algorithm is 
0{n log n + d 3 nt), where n» d and t is the number of cells of the constructed 
triangulation. 

Note that the complexity of BB, without using the method of dynamic de- 
terminants, is bounded by 0(nlogn + d u+1 nt). Recall that t is bounded by 
0(nL<V 2 J) |ZIe95l Sect.8.4], which shows that Alg.[IJ and convex hull algorithms 
in general, do not have polynomial complexity. The schematic description of 
Alg. [T] and its coarse analysis is good enough for our purpose: to illustrate the 
application of dynamic determinant computation to incremental convex hulls 
and to quantify the improvement of our method. See Sect. [4] for a practi- 
cal approach to incremental convex hull algorithms using dynamic determinant 
computations. 

The above results can be extended to improve the complexity of geometric 
algorithms that are based on convex hulls computations, such as algorithms for 
regular or Delaunay triangulations and Voronoi diagrams. It is straightforward 
to apply the above method in orientation predicates appearing in point location 
algorithms. By using Alg. [T] we compute a triangulation and a map of adjoint 
matrices to its cells. Then, the point location predicates can be computed using 
Eq. [4j avoiding the computation of new adjoint matrices. 

Corollary 2 Given a triangulation of a d-dimensional pointset computed by 
Alg. [7| the orientation predicates involved in any point location algorithm can 
be computed in 0(d) time and 0(d 2 t) space, where t is the number of cells of 
the triangulation. 

4 Implementation and Experimental Analysis 

We propose the hashed dynamic determinants scheme and implement it in CH — h 
The design of our implementation is modular, that is, it can be used on top of 
either geometric software providing geometric predicates (such as orientation) or 
algebraic software providing dynamic determinant algorithm implementations. 
The code is publicly available from http://hdch.sourceforge.net 

The hashed dynamic determinants scheme consists of efficient implementa- 
tions of algorithms dynAnv and dyn_adj (Sect.pl) and a hash table, which stores 
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intermediate results (matrices and determinants) based on the methods pre- 
sented in Sect. [3] Every (d — l)-face of a triangulation, i.e., a common facet 
of two neighbor cells (computed by any incremental convex hull package which 
constructs a triangulation of the computed convex hull), is mapped to the in- 
dices of its vertices, which are used as keys. These are mapped to the adjoint 
(or inverse) matrix and the determinant of one of the two adjacent cells. Let us 
illustrate this approach with an example, on which we use the dyn_adj algorithm. 

Example 1 Let A = {a-y = (0,1), a 2 = (1,2), o 3 = (2,1), o 4 = (1,0), a 5 = 
(2,2)} where every point has an index i from 1 to 5. Assume we are in 
some step of an incremental convex hull or point location algorithm and let T = 
{{1,2,4}, {2,3,4}} be the 2-dimensional triangulation of A computed so far. 
The cells of T are written using the indices of the points in A. The hash table 
will store as keys the set of indices of the 2-faces ofT, i.e., {{1, 2}, {2, 4}, {1, 4}} 
mapping to the adjoint and the determinant of the matrix constructed by the 
points oi, a 2 , 04. Similarly, {{2, 3}, {3, 4}, {2, 4}} are mapped to the adjoint ma- 
trix and determinant of 02,03, 04. To insert 05, we compute the determinant 
0/02,03,05, by querying the hash table for {2,3}. Adjoint and determinant of 
the matrix 0/02,03,04 are returned, and we perform an update of the column 
corresponding to point 04, replacing it by 05 by using Eqs. [$| andJ^J 

To implement the hash table, we used the Boost libraries |booj . To reduce 
memory consumption and speed-up look-up time, we sort the lists of indices 
that form the hash keys. We also use the GNU Multiple Precision arithmetic 
library (GMP), the current standard for multiple-precision arithmetic, which 
provides integer and rational types mpz_t and mpq_t, respectively. 

We perform an experimental analysis of the proposed methods. All experi- 
ments ran on an Intel Core i5-2400 3.1GHz, with 6MB L2 cache and 8GB RAM, 
running 64-bit Debian GNU/Linux. We divide our tests in four scenarios, ac- 
cording to the number type involved in computations: (a) rationals where the 
bit-size of both numerator and denominator is 10000, (b) rationals converted 
from doubles, that is, numbers of the form m x 2 P , where m and p are inte- 
gers of bit-size 53 and 11 respectively, (c) integers with bit-size 10000, and (d) 
integers with bit-size 32. However, it is rare to find in practice input coeffi- 
cients of scenarios (a) and (c). Inputs are usually given as 32 or 64-bit numbers. 
These inputs correspond to the coefficients of scenario (b). Scenario (d) is also 
very important, since points with integer coefficients are encountered in many 
combinatorial applications (see Sect.[l}. 

We compare state-of-the-art software for exact computation of the deter- 
minant of a matrix. We consider LU decomposition in CGAL [CGA , op- 
timized LU decomposition in Eigen [GJ + 10| . LinBox asymptotically optimal 
algorithms |DGG + 02] (tested only on integers) and Maple 14 LinearAlgebra 
[Determinant] (the default determinant algorithm). We also implemented two 
division-free algorithms: Bird's |Birll| and Laplace expansion (Poo061 Sect. 4. 2]. 
Finally, we consider our implementations of dynAnv and dyn^adj. 

We test the above implementations in the four coefficient scenarios described 
above. When coefficients are integer, we can use integer exact division algo- 
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rithms, which are faster than quotient-remainder division algorithms. In this 
case, Bird, Laplace and dyn_adj enjoy the advantage of using the number type 
mpz_t while the rest arc using mpq_t. The input matrices are constructed start- 
ing from a random d x d matrix, replacing a randomly selected column with a 
random d vector. We present experimental results of the four input scenarios in 
Tables [TJ [2| [3j [3j We tested a fifth coefficient scenario (rationals of bit-size 32), 
but do not show results here because timings are quite proportional to those 
shown in Table [TJ We stop testing an implementation when it is slow and far 
from being the fastest (denoted with '-' in the Tables). 

On one hand, the experiments show the most efficient determinant algo- 
rithm implementation in the different scenarios described, without considering 
the dynamic algorithms. This is a result of independent interest, and shows the 
efficiency of division-free algorithms in some settings. The simplest determinant 
algorithm, Laplace expansion, proved to be the best in all scenarios, until di- 
mension 4 to 6, depending on the scenario. It has exponential complexity, thus 
it is slow in dimensions higher than 6 but it behaves very well in low dimensions 
because of the small constant of its complexity and the fact that it performs no 
divisions. Bird is the fastest in scenario (c), starting from dimension 7, and in 
scenario (d), in dimensions 7 and 8. It has also a small complexity constant, 
and performing no divisions makes it competitive with decomposition methods 
(which have better complexity) when working with integers. CGAL and Eigen 
implement LU decomposition, but the latter is always around two times faster. 
Eigen is the fastest implementation in scenarios (a) and (b), starting from di- 
mension 5 and 6 respectively, as well as in scenario (d) in dimensions between 
9 and 12. It should be stressed that decomposition methods are the current 
standard to implement determinant computation. Maple is the fastest only in 
scenario (d), starting from dimension 13. In our tests, Linbox is never the best, 
due to the fact that it focuses on higher dimensions. 

On the other hand, experiments show that dyn_adj outperforms all the other 
determinant algorithms in scenarios (b), (c), and (d). On each of these scenarios, 
there is a threshold dimension, starting from which dyn_adj is the most efficient, 
which happens because of its better asymptotic complexity. In scenarios (c) and 
(d), with integer coefficients, division- free performs much better, as expected, 
because integer arithmetic is faster than rational. In general, the sizes of the 
coefficients of the adjoint matrix are bounded. That is, the sizes of the operands 
of the arithmetic operations are bounded. This explains the better performance 
of dyn_adj over the dyn_inv 7 despite its worse arithmetic complexity. 

For the experimental analysis of the behaviour of dynamic determinants 
used in convex hull algorithms (Alg. [I] Sect. |3|, we experiment with four 
state-of-the-art exact convex hull packages. Two of them implement incremen- 
tal convex hull algorithms: triangulation BDH09 implements [CMS93] and 
beneath-and-beyond (bb) in polymake GJOO . The package edd FukQ8] im- 
plements the double description method, and Irs implements the gift-wrapping 
algorithm using reverse search jAviOOj . We propose and implement a variant 
of triangulation, which we will call hdeh, implementing the hashed dynamic 
determinants scheme for dimensions higher than 6 (using Eigen for initial de- 
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d 


Bird 


CGAL 


Eigen 


Laplace 


Maple 


dyri-inv 


dyri-ddj 


3 


16.61 


17.05 


15.02 


11.31 


16.234 


195.38 


191.95 


4 


143.11 


98.15 


71.35 


63.22 


115.782 


746.32 


896.58 


5 


801.26 


371.85 


239.97 


273.27 


570.582 


2065.08 


2795.53 


6 


3199.79 


1086.80 


644.62 


1060.10 


1576.592 


4845.38 


7171.81 


7 


10331.30 


2959.80 


1448.60 


7682.24 


4222.563 







Table 1: Determinant tests, inputs of scenario (a): rationals of bit-size 10000. 
Times in milliseconds, averaged over 1000 tests. Light blue highlights the best 
non-dynamic algorithm while yellow highlights the dynamic algorithm if it is 
the fastest over all. 



d 


Bird 


CGAL 


Eigen 


Laplace 


Maple 


dyri-inv 


dyn_adj 


3 


.013 


.021 


.014 


.008 


.058 


.046 


.023 


4 


.046 


.050 


.033 


.020 


.105 


.108 


.042 


5 


.122 


.110 


.072 


.056 


.288 


.213 


.067 


6 


.268 


.225 


.137 


.141 


.597 


.376 


.102 


7 


.522 


.412 


.243 


.993 


.824 


.613 


.148 


8 


.930 


.710 


.390 




1.176 


.920 


.210 


9 


1.520 


1.140 


.630 




1.732 


1.330 


.310 


10 


2.380 


1.740 


.940 




2.380 


1.830 


.430 


11 




2.510 


1.370 




3.172 


2.480 


.570 


12 




3.570 


2.000 




4.298 


3.260 


.760 


13 




4.960 


2.690 




5.673 


4.190 


1.020 


14 




6.870 


3.660 




7.424 


5.290 


1.360 


15 




9.060 


4.790 




9.312 


6.740 


1.830 



Table 2: Determinant tests, inputs of scenario (b): rationals converted from 
double. Each timing (in milliseconds) corresponds to the average of computing 
10000 (for d < 7) or 1000 (for d > 7) determinants. Highlighting as in Table [lj 
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a 


1 ) • 1 




Eigen 


Laplace 


Linbox 


Maple 


dynAnv 


dyn_adj 


6 


oo 

.23 


3.24 


o c o 

2.58 


1 o 

.15 


1 OO C A 

132.64 


.28 


27.37 


O 1 <7 

2.17 


A 

4 


1.U4 


1 A 

14.01 


i n 

lU.Uo 


.01 


1D4.SU 


l.oD 


( 0. / 


o.oy 


r 

5 


3.40 


45.52 


28.77 


o oo 

2.U2 


367.58 


A CO 

4.52 


176.60 


14.70 


6 


o m 

8.91 


11/1 AC 

114.05 


67.85 


6.16 




/i oo no 

423.08 


OO C G CT 

325.65 


07 A7 

27.97 


7 


20.05 


243.54 


138.80 


42.97 






569.74 


48.49 


o 

8 


40.27 


476.7 4 


OCT O /I 

257.24 








AA /I 01 

904.21 


1 A A 

81.44 


9 


73.90 


815.70 


440.30 








1359.80 


155.70 


10 


129.95 


1358.50 


714.40 








1965.30 


224.10 


11 


208.80 














328.50 


12 


327.80 














465.00 


13 


493.90 














623.80 


14 


721.70 














830.80 


15 


1025.10 














1092.30 


16 


1422.80 














1407.20 


17 


1938.40 














1795.60 


18 


2618.30 














2225.80 


19 


3425.70 














2738.00 


20 


4465.60 














3413.40 



Table 3: Determinant tests, inputs of scenario (c): integers of bit-size 10000. 
Times in milliseconds, averaged over 1000 tests for d < 9 and 100 tests for d > 9. 
Highlighting as in Table [T] 



d 


Bird 


CGAL 


Eigen 


Laplace 


Linbox 


Maple 


dyri-inv 


dyri-adj 


3 


.002 


.021 


.013 


.002 


.872 




045 


.030 


.008 


4 


.012 


.041 


.028 


.005 


1.010 




094 


.058 


.015 


5 


.032 


.080 


.048 


.016 


1.103 




214 


.119 


.023 


6 


.072 


.155 


.092 


.040 


1.232 




602 


.197 


.033 


7 


.138 


.253 


.149 


.277 


1.435 




716 


.322 


.046 


8 


.244 


.439 


.247 




1.626 




791 


.486 


.068 


9 


.408 


.689 


.376 




1.862 




906 


.700 


.085 


10 


.646 


1.031 


.568 




2.160 


1 


014 


.982 


.107 


11 


.956 


1.485 


.800 




10.127 


1 


113 


1.291 


.133 


12 


1.379 


2.091 


1.139 




13.101 


1 


280 


1.731 


.160 


13 


1.957 


2.779 


1.485 






1 


399 


2.078 


.184 


14 


2.603 


3.722 


1.968 






1 


536 


2.676 


.222 


15 


3.485 


4.989 


2.565 






1 


717 


3.318 


.269 


16 


4.682 


6.517 


3.391 






1 


850 


4.136 


.333 



Table 4: Determinant tests, inputs of scenario (d): integers of bit-size 32. Times 
in milliseconds, averaged over 10000 tests. Highlighting as in Table [T] 
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terminant and adjoint or inverse matrix computation) and using Laplace de- 
terminant algorithm for lower dimensions. The main difference between this 
implementation and Alg. [I] of Sect. [3] is that it does not sort the points and, 
before inserting a point, it performs a point location. Thus, we can take ad- 
vantage of our scheme in two places: in the orientation predicates appearing in 
the point location procedure and in the ones that appear in construction of the 
convex hull. We design the input of our experiments parametrized on the num- 
ber type of the coefficients and on the distribution of the points. The number 
type is either rational or integer. From now on, when we refer to rational and 
integer we mean scenario (b) and (d) , respectively. We test three uniform point 
distributions: (i) in the d-cube [—100, 100] d , (ii) in the origin-centered <i-ball of 
radius 100, and (iii) on the surface of that ball. 

We perform an experimental comparison of the four above packages and 
hdeh, with input points from distributions (i)-(iii) with either rational or in- 
teger coefficients. In the case of integer coefficients, we test hdeh using mpq_t 
(hdch_q) or mpz_t (hdch_z). In this case hdch_z is the most efficient with in- 



put from distribution (ii) (Fig. 1(a) distribution (i) is similar to this) while 



in distribution (iii) both hdch_z and hdch_q perform better than all the other 



packages (see Fig. 1(b) I. In the rational coefficients case, hdch_q is competi- 
tive to the fastest package (Fig. [2| . Note that the rest of the packages cannot 
perform arithmetic computations using mpz_t because they are lacking division- 
free determinant algorithms. Moreover, we perform experiments to test the 
improvements of hashed dynamic determinants scheme on triangulation and 
their memory consumption. For input points from distribution (iii) with integer 
coefficients, when dimension ranges from 3 to 8, hdch_q is up to 1.7 times faster 
than triangulation and hdch_z up to 3.5 times faster (Table [4]). It should be 
noted that hdeh is always faster than triangulation. The sole modification of 
the determinant algorithm made it faster than all other implementations in the 
tested scenarios. The other implementations would also benefit from applying 
the same determinant technique. The main disadvantage of hdeh is the amount 
of memory consumed, which allows us to compute up to dimension 8 (Table [4]). 
This drawback can be seen as the price to pay for the obtained speed-up. 

A large class of algorithms that compute the exact volume of a polytope is 
based on triangulation methods ?B EF98j . All the above packages compute the 
volume of the polytope, defined by the input points, as part of the convex hull 
computation. The volume computation takes place at the construction of the 
triangulation during a convex hull computation. The sum of the volumes of 
the cells of the triangulation equals the volume of the polytope. However, the 
volume of the cell is the absolute value of the orientation determinant of the 
points of the cell and these values are computed in the course of the convex hull 
computation. Thus, the computation of the volume consumes no extra time 
besides the convex hull computation time. Therefore, hdeh yields a competitive 
implementation for the exact computation of the volume of a polytope given by 
its vertices (Fig. [I]). 

Finally, we test the efficiency of hashed dynamic determinants scheme on 
the point location problem. Given a pointset, triangulation constructs a data 
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Figure 1: Comparison of convex hull packages for 6-dimensional inputs with 
integer coefficients. Points are uniformly distributed (a) inside a 6-ball and (b) 
on its surface. 




Figure 2: Comparison of convex hull packages for 6-dimensional inputs with 
rational coefficients. Points are uniformly distributed (a) inside a 6-ball and (b) 
on its surface. 



structure that can perform point locations of new points. In addition to that, 
hdeh constructs a hash table for faster orientation computations. We perform 
tests with triangulation and hdeh using input points uniformly distributed 
on the surface of a ball (distribution (hi)) as a preprocessing to build the data 
structures. Then, we perform point locations using points uniformly distributed 
inside a cube (distribution (i)). Experiments show that our method yields a 
speed-up in query time of a factor of 35 and 78 in dimension 8 to 11, respectively, 
using points with integer coefficients (scenario (d)) (see Table El. 
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0.20 
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3 
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35.54 
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4 


0.39 
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0.21 


34.33 


0.82 


35.46 
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4 


0.90 


37.07 


0.47 


35.48 


1.92 


37.17 


260 


5 


2.22 


39.68 


1.08 


38.13 


3.74 


39.56 


500 


5 


5.10 


45.21 


2.51 


43.51 


8.43 


45.34 


260 


6 


14.77 


1531.76 


8.42 


1132.72 


20.01 


55.15 


500 


6 


37.77 


3834.19 


21.49 


2826.77 


51.13 


83.98 


220 


7 


56.19 


6007.08 


32.25 


4494.04 


90.06 


102.34 


320 


7 


swap 


swap 


62.01 


8175.21 


164.83 


185.87 


120 


8 


86.59 


8487.80 


45.12 


6318.14 


151.81 


132.70 


140 


8 


swap 


swap 


72.81 


8749.04 


213.59 


186.19 



Table 5: Comparison of hdch_z, hdch_q and triangulation using points from 
distribution (iii) with integer coefficients; swap means that the machine used 
swap memory. 
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392.55 
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8 
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156.55 


134 


319438 


14.42 


14012.60 
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9 


70 


45.69 


6826 


265874 


0.28 


276.90 


triang 


9 


70 


176.62 


143 


265874 


13.80 


13520.43 
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10 


50 


43.45 


6355 


207190 


0.27 


217.45 


triang 


10 


50 


188.68 
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207190 


14.40 


14453.46 


hdch_z 


11 


39 


38.82 


5964 


148846 


0.18 


189.56 


triang 


11 


39 


181.35 


122 


148846 


14.41 


14828.67 



Table 6: Point location time of IK and 1000K (1K=1000) query points for 
hdch_z and triangulation (triang), using distribution (iii) for preprocessing 
and distribution (i) for queries and integer coefficients. 
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5 Future Work 



It would be interesting to adapt our scheme for gift-wrapping convex hull algo- 
rithms and implement it on top of packages such as jAviOOj . In this direction, 
our scheme should also be adapted to other important geometric algorithms, 
such as Delaunay triangulations. 

In order to overcome the large memory consumption of our method, we 
shall exploit hybrid techniques. That is, to use the dynamic determinant hash- 
ing scheme as long as there is enough memory and subsequently use the best 
available determinant algorithm (Sect.H), or to clean periodically the hash table. 

Another important experimental result would be to investigate the behavior 
of our scheme using filtered computations. 
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